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Abstract. 

Since the mass of the electron is very small relative to 
atomic masses, Thomson scattering of low-energy photons 
(hv <C m e c 2 ) produces thermal Doppler frequency shifts that 
are much larger than atomic Doppler widths. A method is 
developed here to evaluate the electron scattering emissivity 
from a given radiation field which is considerably faster than 
previous methods based on straightforward evaluation of the 
scattering integral. This procedure is implemented in our mul- 
tilevel radiative code (MALI), which now takes full account of 
the effects of noncoherent electron scattering on level popula- 
tions, as well as on the emergent spectrum. Calculations using 
model atmospheres of hot, low-gravity stars display not only 
the expected broad wings of strong emission lines but also ef- 
fects arising from the scattering of photons across continuum 
edges. In extreme cases this leads to significant shifts of the 
ionization equilibrium of helium. 

Key words: radiative transfer - noncoherent Thomson scat- 
tering - non-LTE - stellar atmospheres 



1. Introduction 

The scattering by electrons of low-energy photons (hv <C m e c 2 ) 
is important in astrophysical situations only under the rather 
special circumstance that the electron-scattering opacity is not 
too much smaller than that from all other sources. Thomson 
scattering in these cases has traditionally been regarded as co- 
herent in frequency, and much work on spectrum formation has 
been based on this assumption. However, if the spectrum of the 
object in questions contains features with a width less than or 
comparable to the electron Doppler width, such as strong lines 
and continuum edges, the scattering may have to be treated 
as noncoherent, i.e. the change in the frequency of the scat- 
tered photon by the e lectron Doppler effec t mus t be taken into 



acc ount (Munch |l94g| , Hummer & Mihalas |l967| , Auer & Miha- 
p~968b| ) 



1968a 



1968b|). Although the description of this redistribu- 



las 

tion process as a convolution of the mean intensity J(v) with 
a redistribution function is simple in principle, its inclusion 



in numerical calculations leads to considerable difficulties in 
practice for three reasons. 1) The scale of the redistribution in 
frequency is much larger than the atomic Doppler line widths; 
2) the evaluation of the convolution at each frequency involves 
an integration over a wide band of surrounding frequencies; 3) 
this coupling of frequencies conflicts with the basic strategy of 
approximate lambda iteration (ALI) methods which involve a 
frequency-by-frequency evaluation of the approximate lambda 
operator. 

The transfer equation for two level systems accounting for 
noncoherent electron scattering have been solved in various 
approximation by a number of workers. Early work on this 
problem was based on the form ulation in terms of a "reversing 
layer" by Chandrasekar (1948), who found a solution for the 
case in which the reddening of scattered photons arising from 
the Compton effect was included. Subsequently Munch (1948) 



considered the effects of electron Doppler redistribution, again 
in terms of a "reversing-layer" model, by means of a Fourier 
transform in frequency. This and other early work is summa- 
rized in Chapter 12 of Chandrasekhar fll96r| . Numerical solu- 
tions have been given for li nes in t he atm ospheres of O-type 
stars by Auer and Mihalas ( 
P 
( 



1968a 



1968b), who assumed com- 



plete redistribution for the atomic scattering. Rangarajan et al 
1991) give solutions for parameterized models with both both 



partial and comp lete redistribution i n the line. 

Hillier (1991) and Hamann et al. (1992) have included non- 
coherent electron scattering (NES) in computing the emergent 
line spectra of realistic models of Wolf-Rayet stars, for which 
the effects are clearly observable. However, in these works NES 
is included only during the formal solution for the emergent 
spectra, after the level populations have been fixed from an 
NLTE solution in which NES has not been taken into account. 

In this paper we develop a numerical method for the so- 
lution to radiative transfer problems in which noncoherent 
electron scattering can play an important role. The method 
is based on approximating the electron scattering redistribu- 
tion function by a sum of exponentials. It is shown that only 
two exponential terms give a very accurate approximation. The 
electron scattering emissivity as a function of frequency at each 
point in the atmosphere then can be expressed through the so- 
lution to two simple differential equations in frequency space. 
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The corresponding difference equations can be solved numeri- 
cally by Gaussian elimination. This procedure is easily included 
in any iterative solution of the combined radiative transfer and 
statistical equilibrium equations. 

In previous schemes for evaluating the NES emissivity, the 
computing time for each depth in the atmosphere scales as 
NfNw, where Nf is the number of frequency grid points, and 
Nw is the number of frequency grid points needed to represent 
adequately the width of the electron redistribution function. 
For applications to spectral lines with structure on scales of 
both atomic and electron Doppler widths, Nw can be of order 
of 10-50 or larger. The present method scales more favorably 
as CNf, where C is a small constant, independent of the fre- 
quency grid. It also automatically enforces exact conservation 
of photon number. 

In this paper we describe an implementation of our method 
of treating noncoherent electron scatter ing for the multilevel 

hereafter 



1991 



1992 



ALI code MALI (Rybicki & Hummer 
RHI, RHII). I n con trast to previous methods (Hillier 1991 
Hamann et al. 1992), full account is taken of the effect of NES 



on the level populations, as well as on the emergent spectrum. 
Although the particular application described here is to MALI, 
the method presented here is quite general and should be easily 
adaptable to other codes as well, especially ALI codes. 

It should be made clear that in this latest generalization 
of MALI noncoherent electron scattering is treated not fully 
by ALI, but only by means of ordinary lambda iteration. This 
was not by choice, but resulted from a failure to find suit- 
able approximate electron scattering operators to use in an 
ALI method. In particular, we tried several approximate op- 
erators based on coherent scattering, which is already solved 
non-iteratively in MALI. All of our choices were either unstable 
or had very poor convergence properties. 

Fortunately, for many cases of interest, the mean number of 
photons scatterings due to electron scattering is very moderate, 
of order a few tens. Since the typical number of iterations in an 
ALI solution can be of order many tens to a hundred, the treat- 
ment of NES by ordinary lambda iteration in these cases will 
not substantially change the net number of iterations required 
for a solution. One should also note that numerical accelera- 



tors, such as Ng's ( 1974 ) method, used in MALI, act to improve 
the convergence of ordinary lambda iteration, as well as ALI. 
However, for problems with mean numbers of scatterings of 
order of a hundred or more, the present method is simply not 
suitable, and other methods will have to be developed. 

In Sect. 2 we present the basic description of the electron 
redistribution problem. In Sect. 3 we develop the approximate 
exponential fit to the electron scattering redistribution func- 
tion and show how the emissivity can be found by solving two 
differential equations. We then give results illustrating the ac- 
curacy of our approximation and the basic features of two-level 
transfer problems including electron scattering treated as co- 
herent and incoherent processes. Sect. 4 describes the incor- 
poration of our method into MALI and gives results for a hot 
O-star atmosphere including noncoherent electron scattering. 
As expected, the effects on individual lines were quite notice- 
able, as were the effects within a few electron Doppler widths 
of continua. However, we also found an unexpected and po- 
tentially important effect for certain continua, which, in the 
absence of noncoherent scattering, are strongly in absorption. 
Then scattering of radiation from the stronger continuum be- 
low the edge into the other can cause a substantial anomalous 



ionization of material, and conseqently a substantial change 
in the radiation field across the entire ionizing continuum, not 
just in the neighborhood of the jump. 

2. Basic equations 

The angle-averaged emissivity at frequency v for unpolarized 
electron scattering in the non-relativistic limit can be written 
= (ttE(u), where ctt is the Thomson cross section and 
E(v) is the scattering integral, 



E(y) 



R{u, v)J(u ) dv . 



(1) 



Here J(y) is the mean intensity, and R(y, v') is the unpo- 
larized, angle-averaged electron scattering redistribution func- 
tion. This relation holds at each spatial point in the medium; 
however, we shall suppresss this spatial dependence in the no- 
tation. The noncoherent scattering expressed by Eq. (|l|) is to be 
distinguished from cohere nt sca ttering, for which E(u) — J(y). 

Hummer & Mihalas (1967) derived explicit forms for the 
redistribution functions for non-relativistic electron scattering 
assuming negligible Compton energy shift. It is convenient for 
our purposes to express their results in terms of a modified 
function R(y), 



R(v, u') 



1 



y = 



where the variable y is defined by 
In v — In v 1 



fa 



- i v 



(2) 



(3) 



The quantity (3t is the mean electron thermal speed divided 
by the speed of light c, 



0T = - /^: = i. 8 4xl0- 5 T 1/2 , 
c V m e 



(4) 



where T is the temperature in Kelvin, m e is the electron 
mass, and k is Boltzmann's constant. 
Hummer 



Mihalas (1967) gave two forms for R, depend- 
ing on whether isotropic scattering or the more exact dipole 
scattering is assumed; these are usually distinguished by the 
subscripts A and B, respectively. Their results area: 



Ra(v) = ierfc 



V/4 . 



T erfc T 



(5) 



R B ( y ) = | ierfc iMl - 12 i 3 erfcl|i + 96 i 5 erfci|i 



,11 2 2 1 4 , 1 _ 



y 2 /i 



/3 l 2 . i 4n \y\ r \y\ 

— — y H v )- — - erfc J — - 

K 2 2 y 20 y ' 2 2 



(6) 



The functions i n erfc denote repeated integrals of the error 
function, defined, e.g., in Abramowitz & Stegun ( 1964 ; Sect. 
7.2). These have been expressed in terms of the ordinary error 
function erfc by means of recurrence relations. 



x The occurrences of "erf" in the formulas of Hummer & Miha- 



las (1967) are typographical errors, which should be replaced 



by "erfc" 
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Formulas (0) and (|^), when substituted into Eq. (jij), give 
redistribution fu nction s that differ from those given by Hum- 
mer & Mihalas ( 1967| ), in that the "line center" frequency uq 
does not appear, and the variable y is now defined in terms 
of logarithms of the frequencies v and v' . These changes have 
been made so that the resulting formulas apply for the whole 
spectrum, not just to a region in the neighborhood of a sin- 
gle line. It is appropriate in the non-relativistic limit that the 
redistribution process should depend primarily on the ratio of 
frequencies rather than the differences that appear in the Hum- 
mer and Mihalas formulas. This dependence on frequency ratio 
comes about because the Compton energy shift is negligible in 
the non-relativistic limit and the energy shift of a photon in 
the scattering process is due solely to the Doppler shifts into 
and out o f the rest frame of the electron, (see, e.g., Rybicki & 
Lightman |l979| ; Sect. 7.3). 

It can easily be verified that formulas (H and (|) lead to 
the Hummer & Mihalas forms when applied to line transfer 
for cases of interest in stellar atmospheres, where /3r <C 1, 
implying that the redistribution functions are sharply peaked 
near v ~ v' . An appropriate expansion of y is then 



(7) 



where i/ g can b e taken to be either v or v' , or, as in Hummer 
& Mihalas (1967), the line center frequency vq. Likewise, the 
overall factor in Eq. (^|) can be replaced by 1/(/3t^o)- 

However, for the present work the full logarithmic defini- 
tion Q for y will be used. In fact, it is convenient to introduce 
the logarithmic frequency variable 



(, = log v, 



which will be used instead of the frequency v. Making the 
change of variables in Eqs. (hi) and (□), we obtain 



(9) 



where E(£) and J(£) now denote the scattering integral 
and mean intensity as functions of the variable £. 

Since photon numbers are conserved upon scattering, the 
following integrals should be equal 



E(v)dv f°° J(y)dv 



hv 



hv 



(10) 



which is exactly satisfied for both forms (g) and (|J). 

This normalization condition ( |l2j ) also justifies the use of 
coherent scattering as an approximation when the scale of vari- 
ation of the radiation field is large compared to the electron 
Doppler width. In that case, J(£') can be taken from under the 
integration in Eq. replacing £' by f, and the normalization 
condition implies the coherent result E(£) = J(£). Coherent 
scattering is usually a good approximation for continua, but 
it can fail badly in the neighborhood of lines and continuum 
edges, as we shall see. 

Another quantity of importance is the second moment 



R(y)y 2 dy = 



(13) 



which measures the effective width of the redistribution 
function. This moment has the same value (1/2) for both forms 
(^) and The second moment is particularly important for 
describing noncoherent electron scattering in cases where large 
numbers of scatterings occur, such as in the far wings of lines. 
In this case, the frequency behavior can be well described by 
a random walk that depends solely on the two moments (tl2j) 
and @. 

By contrast, the fine details in the centers of lines depend 
primarily on the absolute value of the function R(y) near y — 0. 
For this reason the value of -R(O) is very important, and also, to 
a lesser extent, its (right) derivative R'(0 + ). These quantities 
differ for the isotropic and dipole cases, and are given by, 



Ra(0) = 

V7T 



£a(0 + ) = --, 



(8) R B (0) = 



(14) 



(15) 



3. The exponential approximation 

An efficient method for the evaluation of the integral in Eq. (^) 
can be based on a simple approximation to the modified redis- 
tribution functions R(y) in the form of a sum of N exponential 
terms: 



1 N 

R(v) = -^2aibiexp(-bi\y\). 



(16) 



We shall call this the exponential approximation, but it 
should be noted that, in terms of true frequencies, the approx- 
imation actually takes the power law form, 



(The intensities here are defined in terms of energy, so one 
must divide by hv to convert to photon numbers.) In terms of 
the logarithmic variable f, this relation may be written, 



1 N 

R(y) = 2 XI aibi 

i=l 



bi I '01 



(17) 



(11) 



where ^< and v > are, respectively, the smaller and larger 
of v and v . 



Substituting Eq. (|9j) and demanding that the resulting 
equation hold for all functions J(£), we find the normalization 
condition 



R(y)dy = l, 



(12) 



3.1. Fitting procedures 

Sets of constants ai and bi can be determined in a number of 
ways, all of which depend on matching some properties of the 
approximate and true functions. In view of the importance of 
conditions (^2|) and (|l3|), we demand that coefficients should 
accurately satisfy the relations 
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E 



N 

E 

i=l 



0, 



(lib; = -. 



(18) 



(19) 



The simplest exponential approximation consists of one 
term (TV = 1), with two independent coefficients ai and &i 
chosen to satisfy both Eqs. ( |ri| ) and (|l^), namely, oi = 1 
and bi — y/2. Since Eqs. ( |l8| ) and ( pj| ) are the same for both 
isotropic and dipole forms for R(y), this approximation does 
not distinguish between them. 

Improved approximations can be obtained by using two 
terms (JV = 2), with four independent coefficients, giving two 
additional degrees of freedom as compared with the single term 
approximation. In view of the importance of the behavior of 
R(y) near y — 0, we would like to use these new degrees of 
freedom to fit the value and slope of the redistribution func- 
tion at the origin, given by Eq. ( |l4[ ) or ([jl]). These conditions 
require that the following two equations be satisfied, 



^flA = 2R(0), 

i=i 

JV 

J2aib 2 i =-2R'(0+), 



(20) 



(21) 



Note that the four conditions expressed by Eqs. ([1 
through (felh are all of the same general form, 



N 

E 

i=l 



(22) 



with different values of m. Approximations for larger N 
constructed by fitting higher order moments and higher order 
derivatives will also satisfy equations of this form. 

Equation ( |2^ ) represents a set of nonlinear equations for 
the unknown values of the coefficients ai andrk. In general, 
numerical methods must be used to solve themcl. We used the 
Newton-Raphson iterative method and confined our search to 
real solutions; complex solutions imply oscillating exponential 
terms that give unacceptable negative values for R(y) for large 

y- 

We attemped to solve the above N — 2 system of equations 
for the isotropic and dipole functions. However, these nonlinear 
equations apparently do not have (real) solutions when both 
the value and slope are given their exact values. Thus we chose 
to fit the value exactly and to fit the slope as closely as possible, 
namely, -0.333 (instead of -0.5) for the isotropic case and -0.6 
(instead of -0.75) for the dipole case. The results obtained in 
this way are summarized in Table 1. 

In Fig. 1 the exact functions R for isotropic and dipole 
case are plotted along with the corresponding exponential ap- 
proximations. The maximum relative errors over the range 



2 In the special case where the values of m are equally spaced, 
equations of the form (B2) can be so lved analytically by Prony's 



Fig. 1. Exact and approximate redistribution functions R(y) for a) 
isotropic and b) dipole 



Table 1. Coefficients for exponential approximation 



Type 


N 


i 


ai 


bi 


A1&B1 


1 


1 


1.000000000000 


1.414213562373 


A2 


2 


1 


5.682790496584 


1.835757192503 






2 


-4.682790496584 


1.986816272751 


B2 


2 


1 


1.690703717290 


1.614249968779 






2 


-0.690703717290 


2.154326524957 



< y < 3 for the N = 1 approximations are about 25% 
for the isotropic and about 17% for the dipole. The maximum 
relative errors for N = 2 approximations over this range are 
both about 10%. For values beyond y — 3 the relative error 
of the approximation increases rapidly, since the true redistri- 
bution function decays asymptotically faster than exponential. 
However, at y = 3 the absolute values of these functions have 
already decreased by about two orders of magnitude from their 
values at y — 0, so the region y > 3 should not be a substantial 
source of error in the emissivities. 

There is reason to expect solutions based on these approx- 
imations to be more accurate than the above maximum er- 
rors might suggest. Note that these approximations have been 
specifically constructed to be exact for the values and moments 
of the redistribution functions that are most critical for the 
solution. Also, the treatment of electron scattering involves 
integration over the redistribution function and this tends to 
average out the errors. 

3.2. Reduction to differential equations 

One of the advantages of the exponential representation is that 
the evaluation of the scattering integral E can be reduced to 
the solution of differential equations. To see this, let us substi- 
tute Eq. (|l|) into Eq. @, leading to 



often produces complex solutions 



method (see, e.g., Hildebrand 1974). However, Prony's method E(£) = ^ ' a%F^ (£), 



(23) 
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where 

1 f°° 

= I exp(-6 I ^ 1 |C-e'l)J(C'X- (24) 

J — oo 

Differentiating twice with respect to £, one easily shows 
that (£) satisfies the second-order differential equation 



d£ 2 



(25) 



One property of this differential approximation is that it 
exactly preserves the normalization condition (p"l|). To show 
this, we integrate Eq. (Eq) over £, which gives 



(26) 



assuming that the first derivative of J 71 ' 1 ' vanishes at both 
endpoints of the region of integration. This will generally be 
true if the frequency range is chosen so that the radiation field 
is sufficiently small at the limits of the range, or by explicit 
choice of boundary conditions (see Appendix A). Now multi- 
plying by at and summing over all i, using Eqs. jl^ ) and (|2^), 
the normalization condition ( [11] ) is recovered. 

The radiation transfer problem can be treated numerically 
by introducing a frequency grid £j, j = 1, . . . , P. All functions 
of frequency are then represented by the set of their values on 
this grid, e.g., Ej — £/(£,■). Eq. (G3) then gives the values of 
electron scattering emissivity on the grid, 



(27) 



in terms of the quantities = F^(£j). These quantities 
can be determined by solving the differential equation ( p5| ) for 
each value of i by means of the well known Feautrier method. 
The details of the numerical method are given in Appendix A. 
A proof that the numerical method also preserves the correct 
normalization ( |ll| ) is given in Appendix B. 

3.3. Relation to the Kompaneets equation 

The differential equation method of the last section bears some 
superficial similarity with the treatment of elec tron s cattering 
using the Ko mpane ets equation (Kompaneets 



and Lightman 



1979 



1957 



Rybicki 



Sect. 7.6), which is closely related to the 
Fokker-Planck equation. However, as we now explain, the two 
methods are very different in their range of validity and in their 
mathematical structure. 

The major advantage of the Kompaneets equation is that 
it includes the major physical effects associated with elec- 
tron scattering: 1) the frequency spreading due to the thermal 
Doppler effect; 2) the frequency shift due to thermal Doppler 
effect (inverse Compton effect); 3) the frequency shift due to 
electron recoil (Compton effect); and 4) stimulated electron 
scattering. Because of this, it is applicable to many aspects 
of Comptonization phenomena. Its major limitation is that it 
can only describe radiation fields that vary slowly on the scale 
of the frequency shift per scattering. In particular, the Kom- 
paneets equation is not appropriate for treating the neighbor- 
hoods of lines or continuum edges. 



The differential equation method of this paper treats only 
the first physical effect above, namely, the thermal Doppler 
spreading effect. However, its major advantage is that it can 
treat radiation fields that vary rapidly on the scale of typical 
frequency shifts. This is essential for treating lines and contin- 
uum edges in typical stellar atmospheres, where, fortunately, 
the other physical effects treated by the Kompaneets equation 
are usually not very important. 

The mathematical distinction between the two methods is 
very striking. Let us compare the Kompaneets equation with 
the N = 1 approximation of this paper, where E(£) = F^'((). 
Then both methods are formulated using second-order differ- 
ential operators. However, the Kompaneets equation expresses 
the emissivity as a differential operator actingon the radiation 
field; in contrast, the differential equation (051) expresses the 
radiation field as a differential operator acting on the emissiv- 
ity. Thus, our second-order operator is not directly comparable 
to the Kompaneets operator but rather to its inverse. It is this 
distinction that allows our method to treat rapidly varying ra- 
diation fields. 

A further mathematical point is that our method can be 
improved by going to larger values of N, with no noticeable 
numerical problems. An analogous improvement of the Kom- 
paneets equation would involve the introduction of higher or- 
der differential operators; however, it is known that this can 
lead to numerical problems, such as instabilities and negative 
solutions. 

3.4- Tests of the method 

In order to test the exponential approximation method, it was 
applied to some simple, parametrized models. A specially writ- 
ten, non-iterative numerical code was developed for this pur- 
pose, so that the results could be evaluated directly, without 
being confused by questions of iterative convergence. We shall 
present one particular such parametrized model here; another 
will be presented in Sect. 4.2. 



Fig. 2. log Fx (arbitrary units) vs. x for parametrized model using 
coherent and various noncoherent scattering approximations 

In Fig. 2 is shown one result of applying the non-iterative 
code. The case is one of pure line scattering in a finite slab of 
total line mean optical depth tl = 2 x 10 4 and total electron 
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scattering optical depth r c = 20. The photon injection rate into 
the line is uniform with depth. The log of the emergent flux F x 
(in arbitrary units) is plotted against x, the frequency relative 
to line center in units of the line Dopper width. The mass of 
the ion is assumed to be that of hydrogen, so that the electron 
Doppler width is ~ 43 in these units. This is an extreme case 
for which line photons are scattered many times (~ 7"e) by the 
electrons, leading to wide wings extending out to many line 
Doppler widths. The two parts of this figure examine the flux 
over very different frequency bands: Fig. 2a shows the details 
near the center of the line, \x\ < 20, while Fig. 2b shows the 
behavior in the far electron-scattering wings, \x\ < 1000. 

Four types of curves are plotted in Fig. 2. The curves 
marked C give the result assuming that the electron scattering 
is coherent. The remaining curves are marked according to the 
type of redistribution function assumed, A for isotropic and B 
for dipole, and also by the number of terms in the exponential 
approximation. 

One sees immediately that the coherent solution is vastly 
different from the other three solutions, which are themselves 
very similar. The behavior for large x shown in Fig. 2b is partic- 
ularly striking, in that all the noncoherent solutions are graph- 
ically indistinguishable in the electron-scattering wings; this 
is explained simply as the result of all three approximations 
having the correct zeroth and second moments of the redis- 
tribution function. The only place where distinctions between 
the three noncoherent cases can be seen is just outside the 
line core, where the differences between the values of i?(0) are 
important. 

The overwhelming impression from Fig. 2 is that the differ- 
ences in the solutions due to different forms of the noncoherent 
approximation are far less important than the differences with 
the coherent approximation. We conclude that any of the vari- 
ous noncoherent approximations will do acceptably well. From 
a theoretical point of view we prefer the dipole form, since it is 
physically more realistic than the isotropic one. Therefore, all 
further numerical calculations reported in this paper will use 
the N — 2 dipole approximation (B2). 



a secondary frequen cy g rid without double points, o n w hich 
the difference Eqs. (Al) with boundary conditions (A3) are 
to be integrated, and the factors weighting J~ and JJ in Eq. 
(A4) are stored. In each iteration, after J„ is evaluated, the 
difference equations are solved at every depth point using the 
modified Feautrier procedure given in Appendix I of RHI and 
the electron-scattering emissivity is evaluated from Eq. (^?|). 

The electron-scattering emissivity and opacity are taken 
together with those of the free-processes as forming the "back- 
ground opacity and emissivity" \c an d Vc which appear in Eqs. 
(2.6) and (2.7) of RHII. These are then updated after every it- 
eration, i.e. are treated by lambda iteration. The cost of this 
operation per iteration scales as the product of the numbers 
of depth and frequency points, i.e. in the same way as inte- 
gration of the monochromatic transfer equations. In practice 
the additional time required for the calculation of the electron- 
scattering emissivity is hardly discernable. However, the addi- 
tional frequency points needed to resolve spectral features on 
the order of the electron-scattering Doppler width do lead to 
moderate increases in the computation time. 

For coherent electron scattering by lambda iteration one 
simply sets Ej = Jj. In RHII, coherent electron scattering was 
treated by direct integration of the monochromatic transfer 
equation in vector form. For a variety of test cases the results 
for coherent scattering calculated from the direct solution and 
by lambda iteration were found to be identical, thus confirming 
the validity of lambda iteration for electron scattering while all 
other processes are included in the ALI procedure. However, 
for the realistic model atmospheres discussed below, the con- 
vergence from LTE initial populations is very slow, and the it- 
erations are started from populations calculated by the direct 
method. The iterative solutions for coherent scattering then 
converges to exactly the same results as the direct solution. 

4-2. Discussion of results 



4. Implementation in MALI 

The form of the approximate lambda iteration method devel- 
oped in RHI and RHII is implemented in the FORTRAN pro- 
gram MALI. Since the publication of RHII the code has been 
generalized to treat line overlaps and an arbitrary number of 
chemical species. 

Although the implementation described here is specific for 
MALI, it should be easily adaptable to other numerical codes, 
especially to ALI codes. 



4-1. Electron scattering by lambda iteration 

The evaluation of the electron-scattering emissivity described 
in Sect. 3 has been implemented in MALI to treat noncoherent 
electron scattering by lambda iteration, while all bound-bound 
and bound-free processes continue to be treated with our ALI 
procedure. 

The frequency grid is constructed with double points at ev- 
ery continuum edge to allow for discontuities in the radiation 
field. On both sides of every edge and line up to twenty addi- 
tional frequency points are included with a spacing of a speci- 
fied fraction of the electron scattering Doppler width at some 
fiducial temperature. A mapping is constructed to and from 



Fig. 3. Electron density and temperature as functions of mass col- 
umn density for Model 2 

To illustrate the type of effects produced by noncoherent elec- 
tron scattering, we have input into MALI photospheric tem- 
perature and density structures from non-LTE models for early 
O-type stars, computed by Dietmar Kunze. These models con- 
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tain only hydrogen and helium; the relative abundances, effec- 
tive temperatures and surface gravities are given in Table 2. 
The model at T e g — 50000 K is represented in Fig. 3; the other 
two are qualitatively similiar. 



Table 2. Model atmosphere parameters 



Model 


T e g 


log g 


N(He)/N(E) 


1 


40000 K 


4.00 


0.100 


2 


50000 K 


3.74 


0.100 


3 


60000 K 


4.20 


0.111 



In our calculations the temperature and total atom density 
as a function of the mass column density were specified and 

held fixed. The original run of electron density was used only 

to calculate the initial LTE populations, and subsequently the 

electron density was calcuated from the ion populations after Fig. 5. Detail from Fig. 4 of the Hen continuum 
every iteration. The converged electron density was found to be 

essentially identical to that of the original model, as were the | 

level populations deep in the atmosphere where the diffusion 
approximation is valid. As distinct from the original models, 
the results obtained here for noncoherent scattering will not be 
in strict radiative equilibrium because of the changes in atomic 
level population. 

We employed simplified atomic models: H° with five bound 
levels, He with 17 and He + with three. Doppler line profiles, 
which were used in the calculation of the atmospheric struc- 
ture, are also assumed here. This is not realistic for the calcula- 
tion of the surface flux, but permits the comparison of coherent 
and noncoherent scattering without introducing further com- 
plications, and should not influence our conclusions. 



Fig. 6. Detail from Fig. 4 of the Hi continuum 



Fig. 4. Surface flux versus frequency for Model 2 computed with 
coherent (dotted) and noncoherent (solid) scattering 

We concentrate now on Model 2, which with T e s = 50000 K 
and log<? = 3.74 lies close to the Eddington limit. The surface 
flux is shown in Fig. 4. The largest differences between the 
results for coherent and noncoherent electron scattering are 

found in the regions around the Hell edge, the Hi edge and Fig 7 _ Dctail from pig 4 of the Lym£m Q lmc 
Lyman a, which are expanded in Figs. 5, 6, and 7, respectively. 



Fig. 8. Ionization fractions as functions of mass column density for 
Model 2 computed with coherent (dotted) and noncoherent (solid) 
scattering 



The gas is essentially completely ionized; the ionization frac- 
tions H°/H, He°/He and He + /He are given in Fig. 8. In Figs. 
4-8 dotted and solid lines represent coherent and noncoherent 
scattering, respectively. 

Most conspicuous in Figs. 4 and 5 are the order of mag- 
nitude increase in the flux in the Hell continuum and the de- 
crease in the flux just below the edge caused by the noncoherent 
electron scattering. These effects arise from the scattering by 
electrons of the intense radiation just below the edge into the 
threshold region where the photoionization cross section has 
its maximum. This increases the ionization rate of He + by an 
order of magnitude, and thus shifts the ionization balance of 
He by the same amount. The resulting decrease in the pop- 
ulations of He and He + are clearly seen in Fig. 8. Note that 
the additional flux in the He II continuum arises not just from 
the photons scattering across the edge, but rather from those 
escaping from the atmosphere because the the He + population 
is reduced. 

The behavior at the H I edge, shown in Fig. 6, is completely 
different. Now the only effect of noncoherent electron scattering 
is appearance of the "spur" just below the edge. This feature is 
caused by scattering of radiation from the Lyman continuum, 
which throughout much of the atmosphere is strongly in emis- 
sion, to the region just below the threshold. This occurs even in 
models for which the surface flux in the Lyman continuum is in 
absorption. As no additional radiation is fed into the threshold 
region of the Lyman continuum, and the amount lost is very 
small, the ionization balance is not significantly affected. 

The narrow emission wings of Lyman a shown in Fig. 7 for 
coherent electron scattering are examples of the Schuster effect. 
Non-coherence reduces this effect by scattering the radiation 
into the very long wings. The slight deepening of the absorption 
feature arises because the radiation trapped in the line core is 
scatterd out into the wings. However, this effect is weak and in 
other cases the core is slightly shallower with non-coherence. 

The general behavior of some of the other spectral lines 
can be summarized as follows. The Lyman j3 line is similiar in 
appearance to Lyman a but less deep. The three He I resonance 
lines and Hell Lyman a are entirely in absorption with broad, 



very shallow wings. Balmer a is weakly in absorption, as are 
most of the He subordinate lines, although a few line of He I are 
weakly in emission with long shallow wings. All emission and 
absorption lines are weaker for noncoherent than for coherent 
electron scattering. 

The above effects are found to a reduced extent in the 
model with T e g = 60000 K, logg = 4.2, and are present only 
weakly in the model with T off = 40000 K, logp = 4.0. The 
ionization equilibrium in the surface layers shifts by roughly 
factors of 3.5 and 1.7, respectively. This shows that the deci- 
sive feature is the essentially complete ionization of the gas, 
so that the weak electron-scattering opacity is not completely 
overshadowed. 

It might also be objected that these phenomena will be 
radically reduced by the lines converging on the series limit. 
Equating the electron Doppler width to the ionization poten- 
tial of a hydrogenic ion shows that the corresponding principal 
quantum number is approximately 



n = 23.3 T. 



-1/4 



(28) 



where T4 = T/(10 4 K). States with n down to at least this 
value will be involved in shifting of photons in or out of the 
continuum. Thus, two or three electron Doppler widths at the 
temperatures considered here corresponds to a state below the 
region of confluence. In other words, radiation present below 
the line merging region will continue to cause additional ion- 
ization. 

4-3. Effect of metal absorption 

The inclusion of metals is likely to reduce both the enhanced 
Hen emission and the Hi spur. However, for hot, low-gravity 
stars near the Eddington limit, absorption by bound-free tran- 
sitions of metal ions should not substantially reduce the shift of 
the ionization balance, as species with appreciable abundances 
have ground state ionization edges lying in the Hell contin- 
uum, and thus will themselves become more highly ionized, 
just as He itself. 

The effect of metal lines, on the other hand, is not easy 
to determine a priori. The dense array of metal lines (arising 
mostly from Fe IV and Fe v) between the H I and H e 11 continua 
illustrated in Fig. 14c of Pauldrach et al (1994) for a plane 
parallel stellar model with T e g = 50500 K and log g = 3.785 
suggests that metal lines could reduce or quench entirely the 
additional ionization of He. More appropriate would be spher- 
ical models includi ng wi nd effects - the so-called unified mod- 
els (Gabler et al ( 1992 ) and references therein). These give 
fluxes in the He 11 continuum larger than thoses in correspond- 
ing plane parallel models by some orders of magnitude, which 
would shift the ionization equilibria of the metals to higher 
stages (the ionization potential of Feiv is only 0.4 eV larger 
than that of Hell). Moreover, deep in atmosphere of a hot 
star, the line blocking shifts to higher frequencies, while the 
scattering of radiation into the Hell continuum persists. The 
combined effects of these mechanisms remains to be investi- 
gated. 

The main effects found here arise from the flow of radiation 
from deep hot parts of the atmosphere in relatively transpar- 
ent parts of the spectrum, followed by scattering by electrons 
into the opaque regions, such as ionizing continuum and strong 
lines. Such effects are not, of course, limited to stellar photo- 
spheres. Winds of hot stars could be affected, as the scattering 
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of radiation in continuua would allow ions to receive ionizing 
radiation from the photosphere which would otherwise be cut 
off by the increasing red-shift of the stellar radiation as seen 
by the material in the wind. Clouds illuminated by an external 
source of ionizing radiation could also experience an increase 
in the degree of ionization by the same mechanism. Whether or 
not these possibilities are realized in any particular case must 
be investigated. The technique developed here should make 
that possible. 
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A. Appendix: The numerical method 

Applying the second-order Feautrier method to Eq. wc 
obtain, 



0T 



F (i) 



2F. 



(i) 



A3-1/2A3 



A j _i/ 2 A j+1/2 



+ 



Aj+i/2Aj 



+if> = Jj, 



for j = 
Aj-1/2 = & 

1 



, (P — 1), where, 



/2 = £j 



A 



(A 



3+1 



i, 



(Al) 



(A2) 



Equation ( ]Al[ ) gives (P — 2) recurrence relations for the P 
unknowns F^ . In order to solve these equations, they must be 
supplemented by appropriate boundary conditions at the ends 
of the grid. A particularly convenient set of boundary con- 
ditions are those of zero derivative, which ensure the proper 
normalization condition (see Eq. [M and the discussion fol- 
lowing it). The zero derivative conditions can be expr essed to 
second-order accuracy using the method of Auer (1967), which 



gives, 



b 2 A 2 
U i "3/2 



+ 1 )F} 



(0 



b 2 A 2 ' 

°i "3/2 



b 2 A 2 p _ 



-F. 



+ 



1/2 



b 2 A 2 p _ 



+ 1)F, 



1/2 



Jl 
Jp 



(A3) 



When applying these difference equations, a problem arises 
at a continuum discontinuity, where the radiation field J has a 
separate left and right limit. In order to represent the radiation 
field near such a discontinuity we assign two values for Jj, 
namely, JJ for the left limit and Jj for the right. In this case, 
the value Jj in Eq. (Al) may be replaced with the quantity 



Aj_ 1/2 J~ +A J+1/2 J+ 



A 



3-1/2 



+ A 



3+1/2 



(A4) 



this procedure is equivalent to assigning separate frequency 
points for the left and right limits and then letting these two 
points approach each other (the proof of this statement will be 
omitted here). 

Equations (|Al| ) and (A3) constitute a tridiagonal system 



of equations for the values of Fj on the grid, which can be 
solved, as usual, by the method of Gaussian elimination. This 
is done for each value of i (i.e., two values for the N=2 expo- 
nential approximations) , and the results summed according to 
Eq. (ji?)) to give the desired values of Ej . The operations count 
for this solution scales linearly with the number of frequen- 
cies P. Furthermore, the coefficient of P is quite small, since 
only simple algebraic operations are involved. This compares 
very favorably with other methods of applying a convolution 
operator to a very unevenly spaced grid. 

The preceding method for computing the emissivity must 
be done separately for each spatial point (depth) in the 



medium. The coefficients of the difference equations (Al) and 
( |A3| ) depend on depth, but solely through the temperature de- 
pendence of (3t, so a good deal of pretabulation of coefficients is 
possible. The computation time necessary to find the emissiv- 
ity at all depths and frequencies scales simply as the number of 
frequency points times the number of depth points. Typically 
these emissivity computations will represent only a moderate 
fraction of the total time needed for the solution of the entire 
radiative transfer problem. 

B. Appendix: Normalization of the method 

Using the discrete equations of Appendix A, it can be shown 
that normalization condition ([y]) holds exactly, providing the 
integrals there are appropriately interpreted in terms of the 
trapezoidal rule. The derivation of this result wi ll in clude the 
possibility of d iscontinuities in J, as given in Eq. (A4). Starting 
with Eq. (Al) (with Jj replacing Jj), we multiply by Aj to 
obtain, 



P W A 



2 J 



2-1/2 ' 



where G% 1/2 = Wr/hfiF^ - Ff>)/A j+1/2 . Likewise, Eqs. 
(|A3j) can be written, 



P (i) A 



2 3 



•3+1/2 



■j-l/2 + ~ Jj A i + i/2 



+ ( G j + l/2 



G 



3-1/2-'' 



(Bl) 



ip«A 3/2 = ±ffA, /a + G« 



4p' i) a 



p-1/2 



= i-F«A 



p-1/2 



G 



(i) 

P-1/2' 



(B2) 



Now summing Eqs. (Bl) for j = 2, . . . , (P — 1), and adding in 
both Eqs. (p2j), we obtain 



E^ i)A ^ + E^w2 

3=2 3=1 
P P-l 

E2 Jj ~ Aj ~ 1/2 + E2 Jj+Aj+1/2 ' 



(B3) 



3=2 



3 = 1 



that is, with an appropriately weighted average of the left 
and right limits of J at the jth point. It can be shown that 



Note that all of the G ^ terms have cancelled out. Rearranging 
summation indices, this can be written, 
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p-1 



E \( F ? + F 3+i)^+i/2 = E \vi + -w^+i/2 (B4) 

3=1 3=1 

Multiplying by a,: and summing over all i, using Eqs. ([L8|) 
and @, we obtain, 

p-i p-i 

£ + B i+1 )A i+1/a = ^ i(J+ + J" +1 )A, +1/2 . (B5) 

3=1 3=1 

This result is the discrete version of the normalization con- 
dition (^), where the integrals have been evaluated using the 
trapezoidal rule over each segment of the grid. 
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